ADCIRC Channel Merging ExampleΒΆ

This notebook demonstrates the complete workflow for generating a mesh with VEWs from channel, land, and background meshes. If you are reading this as an HTML document, it is exported from the jupyter notebook example.ipynb. The shell script example.sh also performs the same operations.

OverviewΒΆ

The workflow consists of 9 main steps:

  1. Combining channel and land meshes with VEWs
  2. Combining channel and land meshes by merging overlapping nodes
  3. Subtracting channel + land coverage from background
  4. Merging channel, land, and subtracted background
  5. Adding land boundaries to the merged mesh
  6. Ensuring elevations of VEW channel nodes to be lower than bank nodes
  7. Adjusting VEW barrier heights to be above the bank nodes
  8. Copying nodal attributes in the background mesh to the new mesh
  9. Updating Manning's n values in the new mesh

Let's go through each step with explanations and visualizations.

Load common modules:

InΒ [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import ScalarMappable
from mpl_toolkits.axes_grid1 import make_axes_locatable
from adcircpy import AdcircMesh
from adcircpy.mesh.fort13 import NodalAttributes
import subprocess

Set common parameters:

InΒ [18]:
vmin = -10
vmax = 10
cmap = None
cmap = 'bwr'
extend = 'both'
lw = 1.0
xlim = [-77.895111, -77.277641]
ylim = [34.022876, 34.557258]

x2yratio = (xlim[1] - xlim[0]) / (ylim[1] - ylim[0])

ylim2 = [34.25, 34.35]
ylimdiff = ylim2[1] - ylim2[0]
xlim2 = [-77.70 - ylimdiff * x2yratio * 0.5, -77.70 + ylimdiff * x2yratio * 0.5]

ylim3 = [34.53, 34.57]
ylimdiff3 = ylim3[1] - ylim3[0]
xlim3 = [-77.3 - ylimdiff3 * x2yratio * 0.5, -77.3 + ylimdiff3 * x2yratio * 0.5]

Define common functions:

InΒ [20]:
def plot_mesh(mesh_file, xlim2=xlim2, ylim2=ylim2, vmin=vmin, vmax=vmax, cmap=cmap, extend=extend, drawvew=False, drawland=False, fulldomain=False, title=None):
    """Plot an ADCIRC mesh with optional VEW and land boundaries at different zoom levels.
    
    """
    
    mesh = AdcircMesh.open(mesh_file)

    if fulldomain:
        figsize = (12.5*3/2-1, 6)
        subplot = [131, 132, 133]
    else:
        figsize = (12.5, 6)
        subplot = [121, 122]
        
    levels = np.arange(vmin, vmax, 0.5)

    cnt = 0

    fig = plt.figure(figsize=figsize)

    # Extent 1
    if fulldomain:
        axes = fig.add_subplot(subplot[cnt])
        axes = mesh.make_plot(axes=axes, vmin=vmin, vmax=vmax, cmap=cmap, extend=extend, levels=levels, title=title)
        axes = mesh.triplot(axes=axes)
        if drawvew:
            axes = mesh.boundaries.weir.gdf.plot(ax=axes, color='lightgreen', linewidth=lw, label='VEW')
            axes.legend()
        if drawland:
            axes = mesh.boundaries.land.gdf.plot(ax=axes, color='brown', linewidth=lw, label='Land')
            axes.legend()        
        axes.set_aspect('equal')
        cnt += 1
            
    # Extent 2
    axes = fig.add_subplot(subplot[cnt])
    axes = mesh.make_plot(axes=axes, vmin=vmin, vmax=vmax, cmap=cmap, extend=extend, levels=levels, title=title)
    axes = mesh.triplot(axes=axes)
    axes.set_xlim(xlim)
    axes.set_ylim(ylim)
    if drawvew:
        axes = mesh.boundaries.weir.gdf.plot(ax=axes, color='lightgreen', linewidth=lw, label='VEW')
        axes.legend()
    if drawland:
        axes = mesh.boundaries.land.gdf.plot(ax=axes, color='brown', linewidth=lw, label='Land')
        axes.legend()
    axes.set_aspect('equal')
    cnt += 1

    # Extent 3
    axes = fig.add_subplot(subplot[cnt])
    axes = mesh.make_plot(axes=axes, vmin=vmin, vmax=vmax, cmap=cmap, extend=extend, levels=levels, title=title)
    axes = mesh.triplot(axes=axes)
    axes.set_xlim(xlim2)
    axes.set_ylim(ylim2)
    if drawvew:
        axes = mesh.boundaries.weir.gdf.plot(ax=axes, color='lightgreen', linewidth=lw, label='VEW')
        axes.legend()
    if drawland:
        axes = mesh.boundaries.land.gdf.plot(ax=axes, color='brown', linewidth=lw, label='Land')
        axes.legend()
    axes.set_aspect('equal')
    cnt += 1

    plt.tight_layout()
    plt.show()
    
    return axes
InΒ [Β ]:
def plot_mannings_n(mesh_file, attrs_file, vmin=vmin, vmax=vmax, cmap=cmap, extend=extend, drawvew=False, drawland=False, fulldomain=False, title=None):
    """Plot Manning's n values from a mesh and fort.13 file.
    
    """
    
    # Load and plot the channel mesh
    mesh = AdcircMesh.open(mesh_file)
    
    attrs = NodalAttributes(mesh)
    attrs.import_fort13(attrs_file)

    get_loc = mesh.nodes.index.get_loc
    tri = [[get_loc(n_id) for n_id in tri] for tri in mesh.triangles.values]

    values = attrs.get_attribute('mannings_n_at_sea_floor')['values'].squeeze()


    if fulldomain:
        figsize = (12.5*3/2-1, 6)
        subplot = [131, 132, 133]
    else:
        figsize = (12.5, 6)
        subplot = [121, 122]
    
    levels = np.arange(vmin, vmax, 0.005)

    cnt = 0

    fig = plt.figure(figsize=figsize)

    # Extent 1
    if fulldomain:
        axes = fig.add_subplot(subplot[cnt])
        axes.tricontourf(mesh.nodes.x, mesh.nodes.y, tri, values, levels=levels, cmap=cmap, extend=extend, vmin=vmin, vmax=vmax)
        axes = mesh.triplot(axes=axes)
        if drawvew:
            axes = mesh.boundaries.weir.gdf.plot(ax=axes, color='lightgreen', linewidth=lw, label='VEW')
            axes.legend()
        if drawland:
            axes = mesh.boundaries.land.gdf.plot(ax=axes, color='brown', linewidth=lw, label='Land')
            axes.legend()        
        axes.set_title(title)
        axes.set_aspect('equal')
        mappable = ScalarMappable(cmap=cmap)
        mappable.set_array([])
        mappable.set_clim(vmin, vmax)
        divider = make_axes_locatable(axes)
        cax = divider.append_axes('bottom', size='2%', pad=0.5)
        cbar = plt.colorbar(mappable, cax=cax, orientation='horizontal')
        cbar.set_ticks([vmin, vmax])
        cbar.set_ticklabels([np.around(vmin, 2), np.around(vmax, 2)])
        cnt += 1
            
    # Extent 2
    axes = fig.add_subplot(subplot[cnt])
    axes.tricontourf(mesh.nodes.x, mesh.nodes.y, tri, values, levels=levels, cmap=cmap, extend=extend, vmin=vmin, vmax=vmax)
    axes = mesh.triplot(axes=axes)
    axes.set_xlim(xlim)
    axes.set_ylim(ylim)
    if drawvew:
        axes = mesh.boundaries.weir.gdf.plot(ax=axes, color='lightgreen', linewidth=lw, label='VEW')
        axes.legend()
    if drawland:
        axes = mesh.boundaries.land.gdf.plot(ax=axes, color='brown', linewidth=lw, label='Land')
        axes.legend()
    axes.set_title(title)
    axes.set_aspect('equal')
    mappable = ScalarMappable(cmap=cmap)
    mappable.set_array([])
    mappable.set_clim(vmin, vmax)
    divider = make_axes_locatable(axes)
    cax = divider.append_axes('bottom', size='2%', pad=0.5)
    cbar = plt.colorbar(mappable, cax=cax, orientation='horizontal')
    cbar.set_ticks([vmin, vmax])
    cbar.set_ticklabels([np.around(vmin, 2), np.around(vmax, 2)])
    cnt += 1

    # Extent 3
    axes = fig.add_subplot(subplot[cnt])
    axes.tricontourf(mesh.nodes.x, mesh.nodes.y, tri, values, levels=levels, cmap=cmap, extend=extend, vmin=vmin, vmax=vmax)
    axes = mesh.triplot(axes=axes)
    axes.set_xlim(xlim2)
    axes.set_ylim(ylim2)
    
    if drawvew:
        axes = mesh.boundaries.weir.gdf.plot(ax=axes, color='lightgreen', linewidth=lw, label='VEW')
        axes.legend()
    if drawland:
        axes = mesh.boundaries.land.gdf.plot(ax=axes, color='brown', linewidth=lw, label='Land')
        axes.legend()
    axes.set_title(title)
    axes.set_aspect('equal')
    mappable = ScalarMappable(cmap=cmap)
    mappable.set_array([])
    mappable.set_clim(vmin, vmax)
    divider = make_axes_locatable(axes)
    cax = divider.append_axes('bottom', size='2%', pad=0.5)
    cbar = plt.colorbar(mappable, cax=cax, orientation='horizontal')
    cbar.set_ticks([vmin, vmax])
    cbar.set_ticklabels([np.around(vmin, 2), np.around(vmax, 2)])
    cnt += 1

    plt.tight_layout()
    plt.show()
    
    return axes

Step 1: Combining Channel and Land Meshes with VEWsΒΆ

First, we combine the channel and land meshes. VEW boundaries are added along the interfaces of the two meshes.

Input mesh (channel):

InΒ [4]:
axes = plot_mesh('input/channel.14', title='Channel Mesh')
No description has been provided for this image

Input mesh (land):

InΒ [5]:
axes = plot_mesh('input/land.14', title='Land Mesh')
No description has been provided for this image

Execute mesh_merger command:

InΒ [6]:
# Combine channel and land meshes with VEWs
# Equivalent to: python -m adcircutils.mesh.mesh_merger input/channel.14 input/land.14 -b vew -o output/ch_la.14 -d 'merged: channel + land with VEWs' 
from adcircutils.mesh import mesh_merger
subprocess.run(['python', '-m', 'adcircutils.mesh.mesh_merger', 
               'input/channel.14', 'input/land.14',
               '-b', 'vew',
               '-o', 'output/ch_la.14',
               '-d', 'merged: channel + land with VEWs'])
Successfully read channel mesh from: input/channel.14
Successfully read land mesh from: input/land.14

Channel mesh info:
Number of nodes: 10505
Number of elements: 17826

Land mesh info:
Number of nodes: 15710
Number of elements: 27617

Starting mesh merge with VEW boundaries...
Creating node mapping...
Combining nodes...
Combined nodes: 26215
Combining elements...
Combined elements: 45443
Finding matching nodes...
Found 2916 matching nodes
Updating boundary node numbers...
Removing matching nodes from boundaries...
Creating final merged mesh...
Finding paired nodes...
Found 2916 paired nodes
Creating VEW boundaries...
Adding VEW boundaries to mesh...
Added 16 VEW boundary segments
Merging endpoints of non-looped VEW boundaries...
Mesh merge complete!

Merged mesh info:
Number of nodes: 26213
Number of elements: 45443

Merged mesh saved to: output/ch_la.14
Out[6]:
CompletedProcess(args=['python', '-m', 'adcircutils.mesh.mesh_merger', 'input/channel.14', 'input/land.14', '-b', 'vew', '-o', 'output/ch_la.14', '-d', 'merged: channel + land with VEWs'], returncode=0)

Resulting mesh with VEWs:

InΒ [7]:
axes = plot_mesh('output/ch_la.14', drawvew=True, title='Channel + Land Mesh with VEWs')
No description has been provided for this image

Step 2: Combining Channel and Land Meshes by Merging Overlapping NodesΒΆ

Next, we combine the channel and land meshes by merging overlapping nodes. This merged mesh is used later for defining the extent of the added mesh for updating Manning's n values.

InΒ [8]:
# Combine channel and land meshes by merging overlapping nodes
# Equivalent to: python -m adcircutils.mesh.mesh_merger input/channel.14 input/land.14 -b merge -o output/ch_la_mgd.14 -d 'merged: channel + land with merged nodes'
subprocess.run(['python', '-m', 'adcircutils.mesh.mesh_merger',
               'input/channel.14', 'input/land.14', 
               '-b', 'merge',
               '-o', 'output/ch_la_mgd.14',
               '-d', 'merged: channel + land with merged nodes'])
Successfully read channel mesh from: input/channel.14
Successfully read land mesh from: input/land.14

Channel mesh info:
Number of nodes: 10505
Number of elements: 17826

Land mesh info:
Number of nodes: 15710
Number of elements: 27617

Starting mesh merge with node merging...
Finding matching nodes...
Found 2916 matching nodes
Creating node mapping...
Combining nodes...
Updating values at matching nodes using channel mesh values...
Adding non-matching nodes...
Combined nodes: 23299
Updating element connectivity...
Combined elements: 45443
Updating boundary node numbers...
Finding boundary edges in final mesh...
Removing non-boundary nodes...
Creating final merged mesh...
Merged 2916 duplicate nodes
Used channel mesh values at 2916 matching nodes
Mesh merge complete!

Merged mesh info:
Number of nodes: 23299
Number of elements: 45443

Merged mesh saved to: output/ch_la_mgd.14
Out[8]:
CompletedProcess(args=['python', '-m', 'adcircutils.mesh.mesh_merger', 'input/channel.14', 'input/land.14', '-b', 'merge', '-o', 'output/ch_la_mgd.14', '-d', 'merged: channel + land with merged nodes'], returncode=0)

Resulting mesh with merged nodes:

InΒ [9]:
axes = plot_mesh('output/ch_la_mgd.14', title='Channel + Land Mesh with Merged Nodes')
No description has been provided for this image

Step 3: Subtracting Channel + Land Coverage from BackgroundΒΆ

Now we subtract the channel and land coverage from the background mesh.

InΒ [10]:
# Subtract channel and land coverage from background mesh
# Equivalent to: python -m adcircutils.mesh.mesh_subtractor input/background.14 output/ch_la_mgd.14 -o output/background_subtracted.14 -d 'subtracted: background - (channel + land)'
subprocess.run(['python', '-m', 'adcircutils.mesh.mesh_subtractor',
               'input/background.14', 'output/ch_la_mgd.14',
               '-o', 'output/background_subtracted.14', 
               '-d', 'subtracted: background - (channel + land)'])
Starting mesh subtraction process...
Finding boundary edges of mesh B...
Found 1153 boundary edges
Splitting boundary into segments...
Found 1 boundary segments
Generating boundary polygon for mesh B...
Boundary polygon generated
Calculating element centroids...
Testing containment of element centroids...
Found 1598517 elements to keep out of 1656634
Creating new mesh with remaining elements...
Renumbering nodes...
Creating final mesh object...
Mesh subtraction complete!
Original mesh: 1656634 elements, 886541 nodes
Result mesh: 1598517 elements, 858058 nodes
Subtracted mesh written to output/background_subtracted.14
Out[10]:
CompletedProcess(args=['python', '-m', 'adcircutils.mesh.mesh_subtractor', 'input/background.14', 'output/ch_la_mgd.14', '-o', 'output/background_subtracted.14', '-d', 'subtracted: background - (channel + land)'], returncode=0)

Resulting mesh:

InΒ [11]:
axes = plot_mesh('output/background_subtracted.14', fulldomain=True, title='Subtracted Background Mesh')
No description has been provided for this image

Step 4: Merging Channel, Land, and Subtracted BackgroundΒΆ

We now merge the channel, land, and subtracted background meshes.

InΒ [12]:
# Merge channel, land, and subtracted background meshes
# Equivalent to: python -m adcircutils.mesh.mesh_merger output/ch_la.14 output/background_subtracted.14 -b merge -o output/ch_la_bg.14 -d 'merged: background_subtracted + channel + land'
subprocess.run(['python', '-m', 'adcircutils.mesh.mesh_merger',
               'output/ch_la.14', 'output/background_subtracted.14',
               '-b', 'merge',
               '-o', 'output/ch_la_bg.14',
               '-d', 'merged: background_subtracted + channel + land'])
Successfully read channel mesh from: output/ch_la.14
Successfully read land mesh from: output/background_subtracted.14

Channel mesh info:
Number of nodes: 26213
Number of elements: 45443

Land mesh info:
Number of nodes: 858058
Number of elements: 1598517

Starting mesh merge with node merging...
Finding matching nodes...
Found 1153 matching nodes
Creating node mapping...
Combining nodes...
Updating values at matching nodes using channel mesh values...
Adding non-matching nodes...
Combined nodes: 883118
Updating element connectivity...
Combined elements: 1643960
Updating boundary node numbers...
Finding boundary edges in final mesh...
Removing non-boundary nodes...
Creating final merged mesh...
Merged 1153 duplicate nodes
Used channel mesh values at 1153 matching nodes
Mesh merge complete!

Merged mesh info:
Number of nodes: 883118
Number of elements: 1643960

Merged mesh saved to: output/ch_la_bg.14
Out[12]:
CompletedProcess(args=['python', '-m', 'adcircutils.mesh.mesh_merger', 'output/ch_la.14', 'output/background_subtracted.14', '-b', 'merge', '-o', 'output/ch_la_bg.14', '-d', 'merged: background_subtracted + channel + land'], returncode=0)

Resulting mesh:

InΒ [13]:
axes = plot_mesh('output/ch_la_bg.14', drawvew=True, fulldomain=True, title='Merged channel, land, and subtracted background mesh')
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/geopandas/plotting.py:482: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  ax.figure.canvas.draw_idle()
/tmp/ipykernel_2363074/121476229.py:66: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  plt.tight_layout()
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/IPython/core/pylabtools.py:170: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)
No description has been provided for this image

Step 5: Adding Land Boundaries to the Merged MeshΒΆ

Next, we add land boundaries to the merged mesh.

InΒ [14]:
# Add land boundaries to the merged mesh
# Equivalent to: python -m adcircutils.mesh.add_land_boundaries output/ch_la_bg.14 -o output/ch_la_bg_lb.14 -d 'merged: background_subtracted + channel + land with land boundaries'
subprocess.run(['python', '-m', 'adcircutils.mesh.add_land_boundaries',
               'output/ch_la_bg.14',
               '-o', 'output/ch_la_bg_lb.14',
               '-d', 'merged: background_subtracted + channel + land with land boundaries'])
Reading input mesh...
Successfully read mesh from: output/ch_la_bg.14

Input mesh info:
Number of nodes: 883118
Number of elements: 1643960
Computing mesh topology...
Found 122572 boundary edges and 122572 boundary nodes

Starting to add land boundaries...
Finding boundary nodes...
Found 122572 boundary nodes
Finding nodes with existing boundary conditions...
Found 120367 nodes with existing boundary conditions
Found 2211 endpoints of unassigned boundary edges
Found 2211 unassigned boundary nodes (including relevant endpoints)
Splitting nodes into continuous segments...
Split into 3 continuous segments
Creating land boundary definitions...
  Created land boundary 1 with 2205 nodes
  Created land boundary 2 with 3 nodes
  Created land boundary 3 with 3 nodes
Adding land boundaries to mesh...
  Appending to existing land boundaries...
Creating new mesh with added boundaries...

Successfully added 3 land boundary segments
Land boundary addition complete!

Output mesh info:
Number of nodes: 883118
Number of elements: 1643960

Saving output mesh...
Successfully saved to: output/ch_la_bg_lb.14
Out[14]:
CompletedProcess(args=['python', '-m', 'adcircutils.mesh.add_land_boundaries', 'output/ch_la_bg.14', '-o', 'output/ch_la_bg_lb.14', '-d', 'merged: background_subtracted + channel + land with land boundaries'], returncode=0)

Resulting mesh:

InΒ [15]:
axes = plot_mesh('output/ch_la_bg_lb.14', drawvew=True, drawland=True, fulldomain=True, title='Merged mesh with land boundaries')
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/geopandas/plotting.py:482: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  ax.figure.canvas.draw_idle()
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/geopandas/plotting.py:482: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  ax.figure.canvas.draw_idle()
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/geopandas/plotting.py:482: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  ax.figure.canvas.draw_idle()
/tmp/ipykernel_2363074/121476229.py:66: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  plt.tight_layout()
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/IPython/core/pylabtools.py:170: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)
No description has been provided for this image

Step 6: Ensuring Elevations of VEW Channel NodesΒΆ

We ensure that the elevations of VEW channel nodes are lower than bank nodes.

InΒ [16]:
# Adjust VEW channel elevations
# Equivalent to: python -m adcircutils.mesh.adjust_vew_channel_elevations output/ch_la_bg_lb.14 -o output/ch_la_bg_lb_adjusted1.14
subprocess.run(['python', '-m', 'adcircutils.mesh.adjust_vew_channel_elevations',
               'output/ch_la_bg_lb.14',
               '-o', 'output/ch_la_bg_lb_adjusted1.14'])
Number of channel nodes altered: 16
Modified mesh saved to: output/ch_la_bg_lb_adjusted1.14
Out[16]:
CompletedProcess(args=['python', '-m', 'adcircutils.mesh.adjust_vew_channel_elevations', 'output/ch_la_bg_lb.14', '-o', 'output/ch_la_bg_lb_adjusted1.14'], returncode=0)

Resulting mesh:

InΒ [21]:
axes = plot_mesh('output/ch_la_bg_lb_adjusted1.14', xlim2=xlim3, ylim2=ylim3, drawvew=True, drawland=True, title='Mesh ensuring channel elevations lower than bank nodes')
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/geopandas/plotting.py:482: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  ax.figure.canvas.draw_idle()
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/geopandas/plotting.py:482: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  ax.figure.canvas.draw_idle()
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/geopandas/plotting.py:482: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  ax.figure.canvas.draw_idle()
/tmp/ipykernel_2363074/2483589787.py:65: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  plt.tight_layout()
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/IPython/core/pylabtools.py:170: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)
No description has been provided for this image

Land boundaries (ibtype=20) should be applied at the ends of un-looped VEWs. This is not necessary when the mesh is used in ADCIRC only simulations. This application of land boundaries at the tips of un-looped VEWs is needed for ADCIRC-SWAN coupled simulations. SWAN gives segmentation fault if the land boundary is not given properly.

Step 7: Adjusting VEW Barrier HeightsΒΆ

We adjust VEW barrier heights to be above the bank nodes.

InΒ [22]:
# Adjust VEW barrier heights
# Equivalent to: python -m adcircutils.mesh.adjust_vew_barrier_heights output/ch_la_bg_lb_adjusted1.14 -o output/ch_la_bg_lb_adjusted2.14

subprocess.run(['python', '-m', 'adcircutils.mesh.adjust_vew_barrier_heights',
               'output/ch_la_bg_lb_adjusted1.14',
               '-o', 'output/ch_la_bg_lb_adjusted2.14'])
Number of barrier heights set: 53464
Modified mesh saved to: output/ch_la_bg_lb_adjusted2.14
Out[22]:
CompletedProcess(args=['python', '-m', 'adcircutils.mesh.adjust_vew_barrier_heights', 'output/ch_la_bg_lb_adjusted1.14', '-o', 'output/ch_la_bg_lb_adjusted2.14'], returncode=0)

Step 8: Copying Nodal AttributesΒΆ

We copy nodal attributes from the background mesh to the new mesh.

InΒ [23]:
# Copy nodal attributes from background mesh to new mesh
# Equivalent to: python -m adcircutils.mesh.attribute_transfer input/background.14 input/background.13 output/ch_la_bg_lb_adjusted2.14 -o output/ch_la_bg_lb_adjusted2.13
subprocess.run(['python', '-m', 'adcircutils.mesh.attribute_transfer',
               'input/background.14',
               'input/background.13', 
               'output/ch_la_bg_lb_adjusted2.14',
               '-o', 'output/ch_la_bg_lb_adjusted2.13'])
/home/sbunya/miniforge/envs/vew1d_3.10/bin/python: No module named adcircutils.mesh.attribute_transfer
Out[23]:
CompletedProcess(args=['python', '-m', 'adcircutils.mesh.attribute_transfer', 'input/background.14', 'input/background.13', 'output/ch_la_bg_lb_adjusted2.14', '-o', 'output/ch_la_bg_lb_adjusted2.13'], returncode=1)

Resulting nodal attribute values (Manning's n):

InΒ [24]:
axes = plot_mannings_n('output/ch_la_bg_lb_adjusted2.14', 'output/ch_la_bg_lb_adjusted2.13', drawvew=True, cmap='turbo', vmin=0.02, vmax=0.18, title='Mannings n values copied from background mesh')
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/geopandas/plotting.py:482: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  ax.figure.canvas.draw_idle()
/tmp/ipykernel_2363074/121476229.py:173: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  plt.tight_layout()
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/IPython/core/pylabtools.py:170: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)
No description has been provided for this image

Step 9: Updating Manning's n ValuesΒΆ

Finally, we update Manning's n values in the new mesh. This step has two substeps:

  1. Selecting channel mesh nodes from the new mesh
  2. Updating Manning's n values at the selected nodes

Step 9.1: Selecting Channel Mesh NodesΒΆ

First, we select the channel mesh nodes from the new mesh.

InΒ [25]:
# Select channel mesh nodes
# Equivalent to: python -m adcircutils.utils.node_selector output/ch_la_bg_lb_adjusted2.14 -o output/ch_la_bg_lb_adjusted2_channel_mesh_nodes.csv -m output/ch_la_mgd.14
subprocess.run(['python', '-m', 'adcircutils.utils.node_selector',
               'output/ch_la_bg_lb_adjusted2.14',
               '-o', 'output/ch_la_bg_lb_adjusted2_channel_mesh_nodes.csv',
               '-m', 'output/ch_la_mgd.14'])
Loading mesh: output/ch_la_bg_lb_adjusted2.14
Mesh loaded: 883118 nodes, 1643960 elements

Selecting nodes by mesh boundary...
Loading boundary mesh: output/ch_la_mgd.14
Generating boundary polygon...
Finding unshared edges...
Found 1153 boundary edges
Creating boundary segments...
Found 1 boundary segments
Creating boundary polygons...
Merging boundary polygons...
Creating buffered boundary...
Finding nodes inside or near the boundary...
Found 26213 nodes inside or near the boundary

Combining selections using union...
Combined selection contains 26213 nodes

Writing 26213 nodes to: output/ch_la_bg_lb_adjusted2_channel_mesh_nodes.csv
Done!
Out[25]:
CompletedProcess(args=['python', '-m', 'adcircutils.utils.node_selector', 'output/ch_la_bg_lb_adjusted2.14', '-o', 'output/ch_la_bg_lb_adjusted2_channel_mesh_nodes.csv', '-m', 'output/ch_la_mgd.14'], returncode=0)

Step 9.2: Updating Manning's n ValuesΒΆ

Now we update Manning's n values at the selected nodes.

InΒ [26]:
# Update Manning's n values at selected nodes
# Equivalent to: python -m adcircutils.utils.node_selector output/ch_la_bg_lb_adjusted2.14 -o output/ch_la_bg_lb_adjusted2_channel_mesh_nodes.csv -m output/ch_la_mgd.14
subprocess.run(['python', '-m', 'adcircutils.nodalattribute.manningsn_extractor',
               'output/ch_la_bg_lb_adjusted2.14',
               'input/ccap_landuse_sample.tif',
               '-s', 'output/ch_la_bg_lb_adjusted2_channel_mesh_nodes.csv',
               'input/ccap_class_to_mn_openwater0.02.csv',
               '-f', 'output/ch_la_bg_lb_adjusted2.13',
               '-o', 'output/ch_la_bg_lb_adjusted2_mn_updated.13',
               '--format', 'fort13'])
Creating nodal extent polygons...
0 / 26213
1000 / 26213
2000 / 26213
3000 / 26213
4000 / 26213
5000 / 26213
6000 / 26213
7000 / 26213
8000 / 26213
9000 / 26213
10000 / 26213
11000 / 26213
12000 / 26213
13000 / 26213
14000 / 26213
15000 / 26213
16000 / 26213
17000 / 26213
18000 / 26213
19000 / 26213
20000 / 26213
21000 / 26213
22000 / 26213
23000 / 26213
24000 / 26213
25000 / 26213
26000 / 26213
Warning: Node 883117 not found in mesh or has invalid neighbors
Warning: Node 883118 not found in mesh or has invalid neighbors
done.
Extracting average mannings n values...
0 / 26200
1000 / 26200
2000 / 26200
3000 / 26200
4000 / 26200
5000 / 26200
6000 / 26200
7000 / 26200
8000 / 26200
9000 / 26200
10000 / 26200
11000 / 26200
12000 / 26200
13000 / 26200
14000 / 26200
15000 / 26200
16000 / 26200
17000 / 26200
18000 / 26200
19000 / 26200
20000 / 26200
21000 / 26200
22000 / 26200
23000 / 26200
24000 / 26200
25000 / 26200
26000 / 26200
done.
Out[26]:
CompletedProcess(args=['python', '-m', 'adcircutils.nodalattribute.manningsn_extractor', 'output/ch_la_bg_lb_adjusted2.14', 'input/ccap_landuse_sample.tif', '-s', 'output/ch_la_bg_lb_adjusted2_channel_mesh_nodes.csv', 'input/ccap_class_to_mn_openwater0.02.csv', '-f', 'output/ch_la_bg_lb_adjusted2.13', '-o', 'output/ch_la_bg_lb_adjusted2_mn_updated.13', '--format', 'fort13'], returncode=0)

Resulting Manning's n values:

InΒ [27]:
axes = plot_mannings_n('output/ch_la_bg_lb_adjusted2.14', 'output/ch_la_bg_lb_adjusted2_mn_updated.13', drawvew=True, cmap='turbo', vmin=0.02, vmax=0.18, title='Mannings n values updated based on land use')
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/geopandas/plotting.py:482: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  ax.figure.canvas.draw_idle()
/tmp/ipykernel_2363074/121476229.py:173: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  plt.tight_layout()
/home/sbunya/miniforge/envs/vew1d_3.10/lib/python3.11/site-packages/IPython/core/pylabtools.py:170: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)
No description has been provided for this image

SummaryΒΆ

In this notebook, we have demonstrated the complete workflow for generating a mesh with VEWs from channel, land, and background meshes. The process included:

  1. Combining channel and land meshes with VEWs
  2. Combining channel and land meshes by merging overlapping nodes
  3. Subtracting channel + land coverage from background
  4. Merging channel, land, and subtracted background
  5. Adding land boundaries to the merged mesh
  6. Ensuring elevations of VEW channel nodes to be lower than bank nodes
  7. Adjusting VEW barrier heights to be above the bank nodes
  8. Copying nodal attributes in the background mesh to the new mesh
  9. Updating Manning's n values in the new mesh
    • Selecting channel mesh nodes
    • Updating Manning's n values at the selected nodes

The final mesh is now ready for use in ADCIRC simulations with proper representation of channels, land areas, and VEW boundaries.